Towards a mesoscopic model of water-like fluids with hydrodynamic interactions 
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We present a mesoscopic lattice model for non-ideal fluid flows with directional interactions, 
mimicking the effects of hydrogen-bonds in water. The model supports a rich and complex structural 
dynamics of the orientational order parameter, and exhibits the formation of disordered domains 
whose size and shape depend on the relative strength of directional order and thermal diffusivity. 
By letting the directional forces carry an inverse density dependence, the model is able to display 
a correlation between ordered domains and low density regions, reflecting the idea of water as a 
denser liquid in the disordered state than in the ordered one. 
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I. INTRODUCTION 



Water is a most common fluid, and yet one still full 
of mysteries. Indeed, water has been reckoned to ex- 
hibit dozens of anomalies, as compared to standard flu- 
ids, primarily the fact of being denser in the liquid than 
solid phase, exhibiting a density maximum at 4°C, i.e. 
above the freezing point (for a vivid non technical de- 
scription, see [? ]). Although a fully comprehensive 
theory of water thermodynamics is still missing, there 
is an increasing consensus that most of these anomalies 
can be traced back to the peculiar nature of the hydro- 
gen bond (HB). The HB interaction plays a vital role on 
structure formation within water. For instance, in wa- 
ter at low temperature, the HB's lead to the formation 
of an open, approximately four-coordinated (tetrahedral) 
structure, in which entropy, internal energy and density 
decrease with decreasing temperature [? ]. The equilib- 
rium thermodynamics, i.e. phase diagram, of water is 
exceedingly rich, and an ab-initio comprehensive analy- 
sis of its properties is beyond computational reach. As a 
result, many models have been developed [? ], including 
lattice ones, which display water-like behavior [???]. 
Such lattice models are typically based on a many-body 
lattice-gas Hamiltonian mimicking the essential features 
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of water interactions, with no claim/aim of/at atomistic 
fidelity [? ]. To the best of our knowledge, these mod- 
els have been employed mostly for the study of equilib- 
rium properties, typically via Monte Carlo simulations. 
Yet, in most phenomena of practical interest, water flows 
and, most importantly, a variety of molecules, say col- 
loids, ions and biopolymers, flow along with it, typically 
in nanoconfined geometries. In the biological context, it 
is well known that the competition between hydrophobic 
and hydrophilic interactions plays a crucial role in affect- 
ing the conformational dynamics of proteins [???]. 
On a larger scale, hydrodynamic interactions are know 
to exert a significant effect on the collective dynamics 
and aggregation phenomena within protein suspensions. 
More generally, hydrodynamic interactions are crucial in 
the presence of confining walls, due to their strong cou- 
pling with resulting inhomogcneities [? ]. 

Based on the above, there is clearly wide scope for a 
minimal model of water behavior, capable of including 
hydrodynamic interactions and geometrical confinement, 
at a mesoscopic level (say tens of nanometers to tens 
of microns). In this respect, a remarkable mesoscopic 
methodology has emerged in the last two decades, in the 
form of minimal versions of the (lattice) Boltzmann ki- 
netic equation [????]. Such lattice Boltzmann 
equations (LBE's) have proven fairly successful in sim- 
ulating a broad variety of complex flows across scales, 
from macroscopic fluid turbulence, all the way down to 
biopolymer translocation in nanopores [? ]. The LB 
approach is mostly valued for its flexibility towards the 
treatment of complex geometries and seamless inclusion 
of complex physical interactions, e.g. flows with phase 
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transitions. Such advantages are only accrued by the 
outstanding computational efficiency of the method, es- 
pecially on parallel computers [? ? ]. To the best of 
our knowledge, however, no LB model for water-like flu- 
ids has been developed as yet. In this paper, we present 
the first preliminary effort to fill this gap. More specifi- 
cally, we develop a new LB framework including HB-like 
interactions, and explore the collective dynamics of the 
mesoscopic water-like fluids, both in free space (homoge- 
neous) and nano-confined (heterogeneous) environments. 

This paper is organized as follows. In Sec.[Il]we present 
the basic elements of our mesoscopic approach, namely: 



in Sec. II A we state the transport equations as described 
by the Lattice Boltzmann model for non-ideal fluids, in 
Sec. Ill Bl we describe the directional interactions mim- 
icking hydrogen-bonds and in Sec. II C we present the 



model for the dynamics of the orientational order param- 
eter. The details of the numerical scheme are provided 
in Sec |III| Our results are presented in Sec. |IV| where we 
first give a qualitative analysis of the mesoscopic model 
(Sec. IV A) and an estimate of the numerical parame- 
ters to be used in the simulations (Sec. IV B). We then 



present numerical results for homogeneous (Sec. IV C ) 
and inhomogeneous (wall-confined) scenarios (Sec. IV D 
and IV E I , without and with hydrodynamic flow. Finally, 
in Sec. |V| we provide an outlook of possible directions of 
future investigation. 



where u = At/r, with the time step At = 1, and r the 
relaxation time towards local equilibrium. The relaxation 
time r fixes the fluid kinematic viscosity v — <? s {t— 1/2), 
where c s is the sound speed of the lattice fluid. Here we 
have taken the mesh spacing Ax = 1 so that c s = 1/ y/3. 

The local equilibrium distribution function is a 
Maxwellian expanded to the second order in the fluid 
velocity [? ? ] and it is described by the distribution 
functions f? q 



i , oCi-U , 9(cj-U) 2 



3U" 

2 c 2 



The weights Wi for the D2Q9 lattice are 



w, 
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i = 



i = 1,2,3,4 
i = 5,6,7,8 



(2) 



(3) 



and we take the propagation speed on the lattice c = 1. 
The macroscopic variables in Eq. ^ are the fluid density 
p, and the fluid velocity U, defined as follows: 



p = p(x,t) =^/i(x,i) 



(4) 



i=0 



and: 



II. MESOSCOPIC MODEL FOR WATER 



A. Transport equations 



U = U(x, t) = u(x, t) + ^^F(x, t) (5) 
p(x,t) 



with 



Our fluid model is based on an extension of the lat- 
tice Boltzmann method for ideal fluids [? ? ]. We use a 
two dimensional (2d) model on the D2Q9 lattice depicted 
in Fig. [T] At each grid node x the velocity distribution 
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FIG. 1: Distribution of the discrete molecular velocities Ci, 
i = 0, ..,8, in the two-dimensional D2Q9 lattice. 

function /i(x, t), i.e. the probability to find a particle at 
location x, moving along the lattice direction defined by 
the discrete speed c,, is evolved according to the kinetic 
equation, with Bhatnagar-Gross-Krook (BGK) approxi- 
mation [? ? ]: 

fi (x + Cl , t + At) - fi (x, t) = -u [fi (x, t) - jf (x, t)] (1) 




p(x, t)u(x, t)=Y^ /*( x > t) c i 



(6) 



i=0 



The forcing term F(x, t) in Eq. |5]) reflects the inter- 
particle interactions. At a microscopic (molecular) level, 
these interactions are given by a combination of van der 
Waals and electrostatic forces, which take into account 
excluded volume, dispersion, directional hydrogen bonds 
and multipolar interactions. In this work, the cohesive 
forces prevailing in the aqueous environment are repre- 
sented by a Shan-Chen pseudo-potential model [? ]. The 
water-water cohesive forces are taken proportional to a 
free parameter Gb, and enter the momentum equations 
^ via the forcing term 



F(x, t) = -G b ip(x, t) ^2 u>ii>(x + Cf, t)c 



where the normalization weights Wi are taken as in 
Eq. ^ and tp is a function of the density, tfj(p) = 
(1 — exp(— p)). Under these conditions, the fluid pressure 
receives a non-ideal contribution from potential energy 
interactions and takes the form: 



P(p) 



(8) 
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This non-ideal equation of state supports a liquid-vapor 
phase transition for Gb < —4 (negative Gb values code 
for attraction) at a critical density p c = In 2 in lattice 
units. Note that hard-core, short-range repulsive interac- 
tions in charge of stabilizing critical phase-separation are 
replaced by a self-driven saturation mechanism, whereby 
the force becomes vanishingly small as p 3> 1, as encoded 
in the exponential dependence of the pseudo-potential ip 
on the density p, and for a uniform density. 



B. Directional interactions 

The main distinctive feature of the present model con- 
sists in the inclusion of directional interactions, aimed 
at mimicking hydrogen bonds at a mesoscopic level. At 
the molecular level, the main effect of hydrogen-bonding 
(HB) between different water molecules is to align a 
donor (hydrogen atom) of a given water molecule to the 
acceptor (oxygen atom) of a neighboring molecule. The 
topology of these links is credited for exerting very pro- 
found effects on the collective behavior of water, and ul- 
timately provides a basis for understanding its numerous 
anomalies. In particular, in three dimensions, at low tem- 
perature HB's tend to favor tetrahedral structures (icy- 
water) which, being less packed than isotropic molecules, 
would explain why ice (the ordered phase) is less dense 
than liquid water (disordered phase). HB's lead to the 
formation of complex and dynamic network structures, 
fueled by the unceasing breaking/formation of new HB's. 
Since these are very short-lived events, we cannot expect 
them to be detectable at a mesoscopic level. At such a 
level, we expect to observe the dynamics of a suitable 
order parameter, which takes a zero value in the disor- 
dered phase and non-zero in the ordered one, the latter 
being promoted and sustained by mesoscopic directional 
interactions. Likewise, at a molecular level, temperature 
takes the form of random noise, promoting jumps be- 
tween the different hydrogen-bond configurations (noise 
breaks the bond). At a mesoscopic level, though, noise 
takes the connotation of a diffusive process, driving the 
system towards a mesoscopically uniform state. 
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FIG. 2: Model of a four-arm molecule. Two arms represents 
donors (d) and two arms are acceptors (a). The vector p = 
(Px,p y ) is oriented in the direction of one donor, as shown in 
the figure. 

To incorporate directional interactions (DI), we en- 
dow the lattice fluid with internal degrees of freedom, 



in charge of responding to orientational forces. These 
degrees of freedom are modeled in terms of four square- 
planar oriented bonding arms, as illustrated in Fig. [2] 
The bonding arms of each "molecule" are either 'donor'- 
like or 'acceptor'-like, to indicate hydrogen or oxygen 
type behavior, respectively. Since the four arms are 
rigid, they are uniquely identified by a single parame- 
ter, e.g. the angle 9 formed by, say, the first arm with 
the x coordinate axis. However, such an angle is not 
a convenient order parameter because of its periodicity. 
As a result, we have introduced a vector order parameter 
P = {Px,Py), pointing in the direction of one donor bond- 
ing arm (see Fig. |2|. The angle 9 relates to the vector 
p through 9 = tan (p y /p x ). The interaction between 
two molecules located, respectively, at x and x,, with 
Y t = x, — x, is described by the following interparticle 
pseudo-potential [? ]: 



VjrB(x,Xi) = 0(p(x))<£(p(xi)) exp 



(r< -Rhb? 



i y^£HB(fc,fcQexp 
■ fc=i k'=i 



nfc • r, 



1 



2-1 

! 1 

2^1 



x exp 



(9) 

where r.; = |rj| and (fii,fc/), with k(k') — 1,..,4, are 
the unit vectors indicating the orientation of the four 
bonding arms of a molecule at grid node x (x,-) (see 
Fig. [3]). Rhb is the equilibrium radial distance of the 
HB's, and o~r and o~e control the radial and angular decay 
of the interactions around such equilibrium (stiffness). 
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FIG. 3: Sketch of two model "molecules" interacting via 
hydrogen-bond potential. 

The matrix e#s(fc, k') is introduced to distinguish be- 
tween donors (arms 1 and 2) and acceptors (arms 3 and 
4). The interaction energy between donors and acceptors 
is equal to the constant £hb < 5 whereas the interaction 
between two donors or two acceptors is set to zero. 
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To include the density-dependent propensity of water 
to form ordered states, which reflects the fact that tetra- 
hedral ordered structures are less compact than isotropic 
disordered ones (the two liquid phases of water), we have 
introduced in Eq. ^ the weight function <j>(p) 



Hp) 



i 



-a(- 



-) 



(10) 



where p max and p m i n are the maximum and minimum 
density of the fluid under the chosen conditions and a > 
is a parameter controlling the range of the density varia- 
tion, which we have set equal to 10 to implement a steeply 
decaying function. 

We can now further justify our choice of the vector or- 
der parameter, p, by noting that the potential in Eq. ([9]), 
with Rhb — V%, has four minima, corresponding to 
9 = 7r/4, 37r/4, 57t/4 and 7tt/A. As a result, cos# alone 
cannot distinguish between these four, leaving two of 
them degenerate, unless sin# is also specified. Clearly, 
the two component vector p does not suffer of such lim- 
itation. 

The magnitude of the total torque on p(x, t) due to the 
HB interaction with its eight neighbors reads as follows: 



C. Dynamics of the order parameter 

The idea underlying the present approach is that the 
mesoscopic description of directional interactions is re- 
flected by the hydrodynamic equation of the vector p. 
Such an equation must include three main effects: macro- 
scopic advection, bond formation due to directional in- 
teractions and bond-breaking due to thermal noise. The 
resulting transport equation reads as follows 



dp 

dt 



u Vp = T + D. 



(16) 



The two terms on the right hand side represent the deter- 
ministic torque, T, due to the hydrogen-bond interaction, 
and thermal diffusion, D, due to translational motion. 
As a result, the water-like fluid is characterized by the 
density p, the velocity U and the rotational vector p. 
The equation ( |16[ ) is evolved concurrently with the LB 
equation ([I]) for the fluid density and velocity. 



Deterministic torque 



A9(*,t) = C T ^c t AtF 2 (x,t) 



(11) 



i=l 



where Ct is a constant to be discussed shortly, q = |c, |, 
At = 1 and Fi (x, t) is the angular component of the force 
between site x and Xj, given by 



F i (x,t) = -^ 



k = l 



udOu 



(12) 



where 9 k = tan 1 (n ky / n kx ) , k = 1, ..,4. Thus 



dV H B(x,Xi) dcosa k 



k = l 



Cidcosoti 



(13) 



with afc the angle in between the unit vector in direction 
fife and the velocity Cj 



COSOtk — 



nfc • Cj 



(14) 



here, n kx — cosO^, n ky = sin6 k . This gives 

(cj — Rhb) 2 



CiFi(x,t) = 0( j o(x))0( j o(x. i ))exp 



4 4 

^ ^ eff B (fc,fc')2(cosafc - l)exp 
' fc=i fc'=i 
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(cosafe - l) 2 



2ai 



■ sm Wfe + Ci y cos Ok 



exp 



2 1 

(15) 



The torque is represented by the force in Eq. (151 which 
is inserted in the equations of motion for p, after pro- 
jection of its components along the x and y directions. 
The two components can be computed as illustrated in 
Fig. |4j The molecule rotates from the initial orientation 
9 = tan^ 1 (p y /p x ) to 9 + AO. The angle A9 is larger 
than for counterclockwise rotation, whereas AO < for 
clockwise rotation. By making use of suitable trigono- 
metric relationships, one obtains: 



T x = -2\p\sin^sin(^ + 9^ 
T y = 2|p|«n4r«w(^r + e )- 



(17) 
(18) 



As discussed previously, the torque is zero when the angle 
is in one of the four degenerate minima k = (2k — 
l)7r/4, k= 1,4. 




FIG. 4: Deterministic rotation, T, of the molecule whose 
initial orientation is 9 = tan" 1 (p y /p x ). Note that |T| = 
2|p|sm(|A#j/2). 
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2. Thermal diffusivity 



The kinetic equation for g x j is: 



At a molecular scale, temperature acts as a HB- 
breaking noise which promotes transitions between dif- 
ferent bond configurations. At a mesoscopic level, such 
bond-breaking effect manifests itself as a diffusion pro- 
cess, driving the system towards an isotropic, disordered 
state (cos 9) = (s'mO) = 0, where brackets stand for en- 
semble averaging over a mesoscopic volume of fluid. As 
a result, in our model thermal diffusion is represented 
in standard laplacian form D = D p Ap, where D p is 
the kinematic diffusivity of the vector p, to be detailed 
shortly. 



3. Hydrodynamic interactions 

Hydrodynamic interactions are automatically included 
by moving the vector p along with the fluid velocity U, 
resulting from the LB advection equation ([!]). These in- 
teractions are expected to play a major role in transport 
phenomena, with the fluid in motion and/or suspended 
bodies moving along with it. In the present study, how- 
ever, we devote special attention to a simpler scenario, 
namely the competition between directional interactions 
and thermal diffusion. Such competition is also analyzed 
in the presence of hydrodynamic flows (Poiseuille) , but 
mostly for illustrative purposes. Quantitative investi- 
gation of the highly complex phenomena resulting from 
the concurrent effects of directional interactions, diffusion 
and hydrodynamic transport, are left to future studies. 



g x ,i(x + Ci,t + 1) - g x ,i(x, t) = 

-w p ( ffa , i (x,t)-^(x,t)) + If 



(21) 



where Tf = WiT x is the x component of the deterministic 
rotation, u p = 1/t p , with t p the relaxation time towards 
local equilibrium, described by the density function 



en 

9 x ,i 



PxWi 



1 



Ci-U 



Vi. 



(22) 



Similar expressions hold for p y by replacing x with y. 

The hydrodynamic limit of this LB model for p yields 
the following continuum macroscopic equations: 



dp 

0t 



u Vp = T + D p Ap, 



(23) 



with Dp = (? s {t p — 1/2). An appealing aspect of LB is 
that one can tune the diffusivity to very small values by 
choosing 



1 



(24) 



with e typically of the order l/N, N being the number 
of lattice sites per linear dimension. 

We note that diffusivity is the emergent manifestation 
of microscopic noise, and acts in such a way as to smear 
out spatial gradients of the vector p. One could still add 
a stochastic source to the rhs of equation (161, in order 
to model random noise. This, however, would not be 
consistent with the mesoscopic aim of the present model. 



III. 



NUMERICAL SCHEME FOR THE ORDER 
PARAMETER 



The equation ( 16 ) for p is essentially an advection- 
diffusion-reaction equation, for which a wide variety of 
numerical methods is available, in particular, the so- 
called, moment propagation method, [? ? ]. For the 
sake of uniformity with the advection equation described 



in Sec. II A and also to secure very low numerical dif- 
fusivity, we have opted for a LB integrator also for the 
equation of the order parameter. Consequently, we inte- 
grate the equation of motion for p = (p x ,p y ) by a lattice 
Boltzmann scheme, applied separately to the two com- 
ponents of the vector. This means that, using the same 
D2Q9 lattice geometry of Sec. |II A[ we define two sets 
of density distribution functions £f s ,i(x, i) and g Vt i(x.,t), 
i = 0, .., 8, respectively for p x and p y , such that 



p x = p x (x, t) = g x ,i (x, t) 

i=0 

8 

p x (x,i)u(x,t) = y^g a ,i(x, t)cj. 



(19) 
(20) 



i=0 



IV. NUMERICAL RESULTS 
A. Analysis of the model 

Before discussing the details of the simulation results, 
a few general comments on the expected qualitative sce- 
nario are in order. The vector p moves with the fluid 
and diffuses at a rate fixed by the diffusivity D p . At the 
same time, it rotates under the effect of the torque asso- 
ciated with directional interactions (DI) with the neigh- 
boring fluid sites. In the absence of DI's, and with the 
fluid at rest u = 0, the order parameter would tend to 
a uniform state (disorder), as dictated by thermal diffu- 
sivity. DI's, on the other hand, tend to place the sys- 
tem on local minima of the interaction potential, thereby 
giving rise to metastable ordered domains. The torque, 
which depends only on angular degrees of freedom, takes 
9(x, y; t) to the local minimum closest to the initial condi- 
tion 9(x, y\t = 0). As a consequence, after an initial tran- 
sient stage, the system settles down into its local minima, 
with no transitions between them. Indeed, transitions be- 
tween different minima can be detected only in the initial 
stage, in which diffusion is still capable of affecting the 
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evolution of 9 towards different minima, because the sys- 
tem is still sufficiently far from equilibrium. Once this 
transient is over, the system freezes into a mctastable 
crystal-like state. In this way, the system attains a state 
of coexistence between short-range uniformity and long- 
range disorder. Within each domain, the system remains 
uniform around the corresponding local minimum angle. 
On the other hand, the spatial distribution of the do- 
mains has a fairly disordered pattern, depending on the 
initial conditions and the (inverse) strength of the diffu- 
sion. A richer scenario, i.e. a longer transient with more 
inter-domain disorder, could be enforced by increasing 
the number of HB arms, so as to enhance the number of 
angular minima, hence their mutual competition. This 
would give rise to a more disordered system, but would 
not change once it falls into the closest minimum, the 
system stays forever. 



B. Numerical estimate of the simulation 
parameters 

The quantitative details of the scenario previously de- 
scribed depend on the specific values of the parameters 
governing the physics of the system, primarily the rela- 
tive strength of Dl's versus diffusion, namely Eh b /kT. In 
the lattice Boltzmann model kT = 1/3. In physical units, 
at ambient temperature T — 300F, kT ~ 2,3KJ/mol. 
For a water dimer, the potential energy of the hydro- 
gen bond is of the order of 60KJ/mol, whereas the free 
energy can be estimated as AF ~ kT/2 [? ]. In the 
present model, the role of the free energy, i.e. the en- 
ergy difference between bounded and unbounded states, 
is played by the parameter ehb in Eq. |9]). Therefore, 
with Ehb = —0.1, we fulfill the condition ehb ~ AF. 

The factor Ct, that, according to Eq. (11), weights 



the effect of the torque, is estimated via the momentum 
equation 



d 2 9 _ d6 



(25) 



where I ~ Ma 2 is the moment of inertia, 7 the drag 
coefficient and T the torque. Thus, upon replacing in- 
finitesimal with finite increments, at steady state, we ob- 
tain A9 ~ jzAt. With M = 1 and a = 1, this yields 



C T oc I/7 (see Eq. (jll])). 

Next, we estimate the relative effect of the torque with 
respect to diffusion. First of all, by Taylor expansion, 

A9 



'111 
tit 



V 



close to the minimum of T, Eq. ( 25 ) gives 
where T' = (dT /d6)\$ e . Therefore, the strength of the 
deterministic rotation can be quantified by the relaxation 
time tt — = 7— ■ On the other hand, the diffusion time 



scale in Eq. ( 16 ) can be computed according to — 
with w the width of the interface, typically of trie order 
of 5 lattice units in LB simulations. The ratio between 



diffusive and deterministic forcing is then given by 
IjD p _ I"/D p a 2 



T'W 2 



(26) 



where the second equality comes by dimensional analysis 
from Eq. ([9]), which leads to T' — \ehb\/o 2 - In our 
numerical simulations, this ratio is changed by tuning 
two parameters, the strength of Dl's via Ehb and the 
diffusivity, via the parameter e. Finally, we set 7 = 1 for 
numerical convenience. 



C. Homogeneous scenario 

We first consider the case of a bulk fluid without 
density-dependent interactions, i.e. where Gb = and 
4>{p) = 1. The heterogeneous case will be considered in 
Sec. |IVD| To study the effect of the different contribu- 
tions to the dynamics of the order parameter, we start 
from the random initial condition shown in the top pan- 
els of Fig. [5j and consider the four different scenarios 
described below. 

Run (a): Dl's only, no diffusion. Diffusion is set to 
zero and the dynamics of the vector p (hence of the 6- 
domains) is dictated solely by the torque. 

Run (b): small diffusivity. Diffusion is turned on 
(e = 0.01), but the torque still dominates the dynam- 
ics [e H b = -0.1). 

Run (c): large diffusivity. Diffusion (e = 0.01) now 
dominates over torque (ehb — —0.0001). 

Run (d): No Dl's, diffusion only. The system dynam- 
ics is entirely dictated by diffusion (e = 0.01). 

The parameters for these simulations are collected in 
Table D 



Run 


e 


Vhb 


SHB 


tt/td 


(a) 





ON 


-1 x 10 _1 


0.0 


(b) 


0.01 


ON 


-1 x 10 _1 


10 -5 


(c) 


0.01 


ON 


-1 x 10" 4 


10" 2 


(d) 


0.01 


OFF 




00 



TABLE I: Simulation parameters. The parameter e fixes 



the numerical diffusivity (see Eq. (24i); tt/td is the rela- 
tive strength of the deterministic torque to the diffusive term 
(see Eq. ( 26 1 ) . The drag coefficient is 7 = 1, with Ct = 1/7 
in En. (111. When the HB potential Vhb is switched on (see 
Eq. 9k, we used the values ag = or = 0.12, R HB = \/2 and 
we set 4>(p) = 1 to first consider the homogeneous scenario. 
The value of ehb is given in the table. The lattice sizes were 
N x = N y = 64 in all cases. 



The equilibrium configurations for 8 in Run (a) and 
Run (b) are shown in Fig. [5J together with the (ran- 
dom) initial configuration for all runs presented in this 
section. As expected, the directional interactions drive 
the system towards a disordered collection of uniform, 
ordered domains, i.e. within each domain 9 takes the 
same value, but domains with different 9 coexist within 
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FIG. 5: Configuration of 6 in a 64 2 homogeneous system, 
shown as contour plot (left), detail of the orientation of the 
four-armed "molecules" (middle, only donor arms are repre- 
sented), and histogram (right). Starting from the same initial 
(n=0) random configuration shown in the top panels, we re- 
port the equilibrium configurations in the following two cases; 
Central panels: Run (a), the scenario with e HB — —1 x 10" 1 
and no diffusion in the dynamics of p (i.e. e = 0). Equi- 
librium is reached after n w 10 3 steps. Bottom panels: Run 
(b), the scenario with €hb = —1 X 1CP 1 and diffusion with 
e = 0.01. Equilibrium is reached after n w 10 4 steps. Note 
that in both cases the histograms are peaked around the four 
equilibrium values = 7r/4,37r/4,57r/4, 7ir/4, but the size of 
the domains increases when diffusion is turned on. Also, the 
peaks of the histogram broaden with increasing diffusion, in- 
dicating thicker domain boundaries. 



a disorderd global pattern, in which all four minima are 
roughly equally populated (as shown by the histograms). 
As discussed previously, the spatial distribution of the 
ordered domains reflects the initial condition, since 9 is 
attracted to the minimum closest to its initial value, and 
once there, it stops evolving because the system is deter- 
ministic (no stochastic fluctuations). This corresponds 
to an arrested-coarsening scenario, as observed in many 
slow-relaxing materials, including water. In the absence 
of diffusion, coarsening is virtually quenched, and the 
resulting domains are very small and very close to the 
initial condition. 

When diffusion is turned on, the system experiences 
additional freedom to evolve away from the initial condi- 
tion. As a result of diffusive transport, coarsening can 
now take place, as indicated by the increased size of 
the uniform domains. Since DI's are still dominant, the 
global pattern remains disordered. As expected, the size 
of the ordered domains depends on the relative strength 
of the HB to the thermal diffusion. Whenever diffusion 
dominates over the HB formation, coarsening can pro- 
ceed up to the point where all domains merge into a single 



one, i.e. a uniform value of 9 all over the region occupied 
by the fluid. This behavior is illustrated in Fig. [6j The 
final value of 9 attained by the fluid corresponds to one 
of the four equilibrium angles of Vhb, and its specific 
value is selected by the initial condition. Also note that 
during its relaxation to global equilibrium, the system 
first reaches local equilibrium (shown by the peaks in the 
histogram around the metastable values of 9), and then 
the ordered domains grow and merge, with the largest 
one absorbing the others. Although coarsening can pro- 
ceed up to complete uniformity, the system still keeps 
track of HB interactions, in that the final 9 cannot just 
take any value, but only one of the four degenerate min- 
ima. This means that, even upon averaging over ini- 
tial conditions, the system remains non-ergodic, meaning 
by this that the probability distribution function (pdf) 
P{9) is given by a combination of the four Dirac's deltas, 

p {0) = | Et=i s ( e - k), with e k = {2k - 1) % k = 1,4. 

This scenario lies at other extreme as compared to the 
no-DI's situation. As illustrated in Fig.[7j diffusion drives 
the system to a uniform state, with a relaxation dynam- 
ics similar to the case with weak directional interactions 
(cfr. with Fig. [6]). However, in this case, the interme- 
diate and final values of 9 can take any value between 
[0, 2w], depending on the initial conditions. On a statis- 
tical basis, this is very different from the case of weak 
DI's because ensemble averaging over initial conditions 
would now produce a uniform pdf, P{9) = 

We wish to emphasize that in the case of diffusion- 
dominated scenarios, the steady-state solution of Eq.(23) 
reads simply as p x = p cos 9 — const and p y = p sin 9 = 
const, where const means constant in time and space. 
This means that both the magnitude p and the orien- 
tation 9, are uniform in space. Such uniformity does 
not correspond to any microscopic "symmetry-breaking" , 
but simply reflects the diffusive structure of the meso- 
scopic equation of motion of the parameter p. Indeed, 
in the diffusion-dominated scenario, isotropy is recovered 
by averaging upon initial conditions. 

The results discussed in this section show that the 
present LB model with directional interactions is capable 
of sustaining long-lived, metastable states in the form of a 
disordered collection of uniform domains, whose asymp- 
totic size and spatial distribution is controlled by the rel- 
ative strength of directional interactions versus diffusion. 



D. Heterogeneous scenario: hydrophobic effects 

The above results pertain to a homogeneous scenario, 
with no hydrodynamic motion and away from solid 
boundaries. However, most important water-mediated 
phenomena take place in the proximity of solid bound- 
aries, where heterogeneity plays a major role [? ? ]. 

In view of future applications, it is therefore impor- 
tant to include these effects in our model. To this pur- 
pose, we have leveraged the LB capability to include non- 
ideal interactions (hydrophobic/hydrophilic) through 




FIG. 6: Run (c). Snapshots of 9 (left panels) and corre- 
sponding histograms (right panels) at different steps n dur- 
ing the simulation for a 64 homogenous system with e HB = 
— 1 x 10 -4 . At equilibrium the histogram is peaked at the 
value 8 = 7n/4, which corresponds to one of the four minima 
of Vhb- The value of the attained minimum depends on the 
initial condition (here, same as in Fig. [5]). 



density-dependent pseudo-potentials (Shan-Chen model, 
Eq. Q). In particular, by adjusting the strength, Gb, 
of attractive fluid-fluid interactions in the bulk, we can 
generate an effective repulsion at walls, since a near- wall 
fluid layer would experience an attractive force from the 
second inner layer in the fluid and no force (or a smaller 
one) from the solid layer at the wall, thus resulting in 
an effective repulsion from the wall (hydrophobic effect). 
This strategy has been used by many authors in LB sim- 



FIG. 7: Run (d). Snapshots of 6 (left panels) and corre- 
sponding histograms (right panels) at different steps n during 
the simulation of a homogenous system of 64 2 without DI's. 
The final value of 6 is random, and it depends on the initial 
condition (here, same as in Fig. [HI). 



ulation of confined microfluids, and shown to give rise 
to density depletion layers near the wall (see [? ] and 
references therein). Coupling the effect of these deple- 
tion layers with the density-dependent prefactor <f){p) of 
Eq. (10), permits to modulate the fluid structure as a 



function of distance from the solid walls. 

For our tests, we have chosen a channel of size N x = 32 
and N y — 128 with two layers of solid sites perpendicular 
to the x-direction, located at x = 2.5 and x = N x — 1.5 
(for computational reasons, we use two buffer layers and 
effective walls lie in between the second buffer and the 
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first fluid layer). We used the half-way bounce back 
scheme for no-slip boundary conditions, i.e particles ex- 
iting the fluid domain are bounced backed into the fluid 
with the opposite velocity, along both normal and tan- 
gential directions. We set the parameter for the Shan- 
Chen interaction (see Eq. (JtJ) ) to Gb = —3.8, a value 
slightly below the critical threshold leading to liquid- 
vapor phase separation. We used the value Ehb = — 0.1 
in Vhb- The initial condition was set to p(x) = In 2 
everywhere within the channel. The density at the 
wall sites was fixed to p wa \\ = 0.55 < In 2 and we set 
Pmin = Pwall- With these parameters, the maximum den- 
sity (reached at the center of the channel, as illustrated 
in Fig. [8]), was p max = 0.72, corresponding to a density 
depletion ratio of about 25%. We first consider the static 
scenario (no net flow) with tt/tjj — 10 -6 and e — 0.001. 

The averaged density profile p(x) — Yly=i p( x > v) 
as a function of the distance from the wall is shown in 
Fig. [8] (top). From this fi gure, a depletion layer extending 
about S — 4Ax (4 lattice sites) away from the wall is well 
visible. The width of this layer is of the same order of 
magnitude of the interaction range, namely y/2 Ax in the 
present nine speed lattice model. 

To inspect the ordering effect, we have computed the 
correlation function of the 8 domains, namely 



0.75 



C e (x;r y ) 



< 9(x,y + r y )9(x,y) > 
< 9(x,y) 2 > 



(27) 



where the brackets stand for ensemble average over initial 
conditions and along the y-axis. The correlation func- 
tion as a function of r y , at various locations away from 
the wall, x = 4,8,10,16 (in lattice units), is shown in 
Fig. 8] (bottom). Near the wall, correlations decay ba- 
sical y within one lattice site, because the strong effect 
of the torque dominates over diffusion effects, quench- 
ing the system around the minima of Vhb closest to the 
initial value of 9. On the other hand, at the center of 
the channel, correlations persist over larger distances, 
because diffusive effects lead to ordering within larger 
domains. The contour plot of the ^-domains is shown 
in Fig. [9] (left panels) , together with the corresponding 
histogram. 



E. Flowing systems 

In the Introduction, we have emphasized the crucial 
role played by hydrodynamic interactions in water trans- 
port phenomena under confinement [? ]. We have also 
highlighted the remarkable capability of LB methods to 
include hydrodynamic effects in the study of the dynam- 
ics of complex flows, including polar ones [? ? ], at virtu- 
ally no extra computational cost. In this section, we pro- 
vide just a qualitative example of such capability, leaving 
an in-depth analysis and applications to future publica- 
tions. We consider the heterogeneous scenario discussed 
in the previous section, in the presence of a Poiseuille flow 
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FIG. 8: Density profile (top) as a function of the distance from 
the walls located at x — 2.5 and x = N x — 1.5 and correlation 
function of 9 (bottom) as a function of r y at various locations 
x, from close to the wall to the center of the channel. 



between the two confining walls. Specifically, u x — and 
u y (x) — AUo(x/W)(l — x/W), where Uq is the centerline 
speed and W the channel width. The relevant dimcn- 
sionless parameter, measuring the strength of convection 
versus diffusion, is the Peclet number Pe = UoW/D p . 
With e = 0.001, hence D p = 0.001/3, W = 32 and 
Uq = 0.02, we have Pe ~ 2000, indicating a strong dom- 
inance of advective effects over diffusive ones. However, 
advection still is three orders of magnitude slower than DI 
interactions, which therefore remain the dominant mech- 
anism. Looking at the typical timescales, we have that 
tt/td = 10~ 6 and we can estimate the ratio T a( i v / tjj as 
~ 1/Pe. Note that the typical advection time r a( j„ takes 
an intermediate value between the torque and diffusion 
times. 

The main qualitative effect of advective motion, and 
particularly of the shear s xy = du y /dx ~ Uq/W, is 
to promote mechanical deformation of the 9 domains, 
as one can appreciate from Fig. [9] (right panels). This 
deformation is accompanied by a mild extra-mixing, as 
witnessed by the broadening of the 9 histograms in the 
bottom panel of the same figure. This is in line with 
the qualitative idea of advection as an additional mixing 
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FIG. 9: Left panels: Contour plot and corresponding his- 
togram of the 9 domains in a heterogeneous system. The 
histogram is nearly equi-distributed among the four degen- 
erate minima of the HB potential. Right panels: Contour 
plot and corresponding histogram of the 9 domains in a het- 
erogeneous system with flow advection. The histogram is still 
nearly equally distributed among the four degenerate minima, 
but with an enhanced spreading due to advective motion. 

mechanism, which, in the limit of high Peclet numbers, 
becomes correspondingly more effective than diffusion in 
competing with the ordering effect of directional interac- 
tions. 

As stated above, a quantitative inspection of these 
complex phenomena will be left for future studies. Nev- 
ertheless, Fig. [9] conveys an idea of the kind of complexity 
that can be tackled by the water-like LB model presented 
in this work. 



V. CONCLUSION AND OUTLOOK 

Summarizing, we have developed a mesoscopic lattice 
Boltzmann model for water-like fluids. By water-like, 



we imply that the fluid includes directional interactions 
(DFs) mimicking (some) generic features of hydrogen- 
bonds, particularly the tendency to form ordered do- 
mains, with selected bond angles, against the action of 
thermal noise. The competition between the ordering ef- 
fect of directional interactions against thermal disorder 
is shown to lead to the formation of complex patterns of 
the orientational order parameter, consisting of a disor- 
dered collection of ordered domains (i.e., in each domain 
the order parameter takes a uniform value). In addition, 
we have included non-ideal interactions which permit to 
realize heterogeneous density depletion layers near solid 
walls. By letting the DI's carry an inverse density depen- 
dence, the model is able to incorporate a built-in corre- 
lation between ordered domains and low density regions, 
reflecting the idea of water as a denser liquid in the dis- 
ordered state than in the ordered one. 

The extension of our model to the three dimensional 
(3d) case does not present any conceptual barrier. The 
3d water molecule would be represented by four bonding 
arms, arranged in a tetrahedral symmetry [? ]. This 
geometry requires two independent orientation angles for 
each lattice site. As far as the lattice is concerned, a 
possibility to accommodate the tetrahedral geometry is 
the standard cubic 15-speed lattice commonly used in 3d 
lattice Boltzmann simulations. 

This paper sets a methodological starting basis for LB 
models of water-like fluids. A fully quantitative valida- 
tion against microscopic models remains as a task for 
the future, warranting a separate investigation on its 
own. The model is expected to prove mostly valuable 
for the study of nanoscopic water flows with suspended 
molecules and/or charged ions, both in standalone [? ] 
and/or multiscale scenarios [? ? ]. 
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